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ABSTRACT 



Aims. In view of the substantial uncertainties regarding the possible dynamics of the dark energy, we aim at constraining the expansion 
rate of the universe without reference to a specific Friedmann model and its parameters. 

Methods. We show that cosmological observables integrating over the cosmic expansion rate can be converted into a Volterra integral 
equation which is known to have a unique solution in terms of a Neumann series. Expanding observables such as the luminosity 
distances to type-la supernovae into a series of orthonormal functions, the integral equation can be solved and the cosmic expansion 
rate recovered within the limits allowed by the accuracy of the data. 

Results. We demonstrate the performance of the method applying it to synthetic data sets of increasing complexity, and to the first- 
year SNLS data. In particular, we show that the method is capable of reproducing a hypothetical expansion function containing a 
sudden transition. 
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1. Introduction 

Astronomical measurements can constrain cosmology through 
two functions, the cosmic expansion rate and the growth rate 
of cosmic structures. If space-time is on average homogeneous 
and isotropic and topologically simply connected, it must be de- 
scribed by a Robertson- Walker metric characterised by a time- 
dependent scale factor a{t). General relativity enters when the 
scale factor is to be related to the energy content of the uni- 
verse. However, the geometry of space-time, in particular the 
distance measures, are determined already once the scale factor 
and its first derivative are specified by the expansion function 
H(a) := a I a. 

The dynamics of structure growth requires more than that. 
Linear structure growth is commonly based on the Newtonian 
approximation to the continuity and Euler equations, together 
with Poisson's equation. They provide a valid foundation if 
Minkowskian space-time is a good local approximation on the 
scales of the structures considered. Then, the total inhomoge- 
neous energy density drives its further evolution in time, and the 
relation of time to observables such as redshift again requires the 
expansion function H(a). 

We are writing this to emphasise that the expansion function 
is the central mathematical object underlying all cosmological 
constraints, augmented by the assumption of local Newtonian 
dynamics if structure growth in the late universe is to be in- 
cluded. This suggests that measurements of the expansion func- 
tion itself, without any reference to Friedmann models, should 
be possible and much more fundamental than the common con- 
straints of cosmological parameters entering the expansion func- 
tion once the density contributions to the Friedmann models are 
specified. 
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This approach may become interesting for instance in the 
context of dark energy, for which our lack of understanding al- 
lows an only poorly constrained, vast range of possible mod- 
els. Usually, dark-energy models are characterised by a small 
set of parameters and placed into the cosmic expansion rate 
through Friedmann's equation, substituting the conventional 
cosmological-constant term. While this is doubtlessly reason- 
able when testing specific dark-energy models, the question re- 
mains interesting what can be inferred on the cosmic expansion 
rate from observations without any reference to a specific model 
for the energy content of the universe and how it may affect cos- 
mic dynamics. 

The importance of a model-independent reconstruction of 
the cosmic expansion rate from luminosit y distance d a ta has 
been largely discussed in the literature. In Starobinsky (1998) 
the relations between the observational data and the expansion 
rate are presented, and several different techniques have been de- 
veloped since then to appropriately tre at the data in or der to per- 
form s uch a reconstruc ti on (see, e.g., Huterer & Turner ( 1999, 
2000); iTegmarkl (120021) : I Wang & Tegmarkl d2005l) ) ."Recent re- 
co nstruction techniques with appli cations to data can be found 
in iDalv & Diorgovskil (|2003, 2004), where luminosity-distance 
data from type-la supernovae are combined with angular- 
diame ter distances inferred from radio galaxies, iFav & Tavakoll 
(2006), where constraints fro m baryon acoustic osci ll ations are 
added to supernova data, and lShafieloo et alj (120061) : IShafielool 
(2007). All these methods employ a smoothing procedure in red- 
shift bins. Also principal component analysis has been used to 
recon struct the dark energy equat i on of state as a function of red- 
shift dHuterer & Starkmanl 120031: ICrittenden & Pogosianl 120051: 
ISimpson & Bridlell2006l) 

We are here proposing a method aiming at a direct deter- 
mination of the cosmic expansion rate without reference to any 
specific model for the energy content of the universe. As pro- 
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posed and used here, it only employs the minimal assumptions 
(1) that the universe is topologically simply connected, homo- 
geneous and isotropic on average, and (2) that the expansion 
rate is reasonably smooth. We introduce the method in Sect. 2 
and demonstrate its performance with synthetic and real type- 
Ia-supernova data in Sect. 3. We finish with the discussion in 
Sect. 4. 

2. Method 

2.1. Preliminary remarks 

The cosmological standard model, which was defined and tightly 
constrained approximately during the last decade, is based on 
the combination of a multitude of cosmological measurements. 
Perhaps unduly simplifying the picture, the typical angular size 
of CMB temperature fluctuations constrains the overall spa- 
tial curvature and thus the sum of all energy density contribu- 
tions (Spergel et al. 2007), the large-scale galaxy power spec- 
trum (Tegmark et al. 2004) and the evolution of galaxy clusters 
dAllen et al.l|200 3, 2004) constrain the total matter density, pri- 
mordial nucleosynthesis constrains the total density of baryonic 
matter (Kneller & Steigman 2004), and difference between the 
total energy density and the matter density is attributed to the 
cosmological constant or the dark energy. 

Type-la supernovae stand out in this context because they 
allow in principle to directly measure the accelerated cos- 
mic expansion, which in the context of Friedmann models re- 
quires a dominant energy contribution with negative pressure. 
Measurement accuracies now seem high enough to not only 
favour accelerated expansion, but to exclude decelerated expan- 
sion. Besides, type-la supernovae constrain cosmology purely 
geometrically because they allow measuring how the luminosity 
distance depends on red shift or, equivalentl y , on the scale facto r 
a. For recent result s seelAstier et a l. (2006) JRie"ss et all (l2007l) : 
IWood-Vasev et all (2007); Davis et al. (2007). 

Without loss of generality, we write the expansion function 
H(a) in the form H(a) - HoE(a), introducing an expansion 
function E(a) arbitrarily normalised to unity at present, a — 1 . 
Already the Robertson- Walker metric, before its specialisation 
to a Friedmann metric, shows that the angular-diameter distance 
is given by 



D A (a) = af K [x(a)] , 

with the comoving angular-diameter distance 

(sin* (K=l) 

f K (x) = \X (K = 0) 

[ sinh^ (K = -1) 

and the comoving distance 



X( a ) = 



H 



f 1 da' 

J a a' 2 E(a') 



(1) 



(2) 



(3) 



Due to Etherington's relation, which holds for any space-time, 
the luminosity distance is 



D L (a) 



f K [x(a)l 



(4) 



In the common approach to constraining cosmological pa- 
rameters and specifically the dark energy, a particular Friedmann 
model is adopted, whose expansion rate is written as 



E(a) 



— + — 3 — + iiqO^(fl) 



1/2 



(5) 



in terms of the present-day density parameters Q( r , mj k, q )0 for 
the radiation, matter, curvature and dark-energy contributions, 



where Qko = 1 — ^rO - ^mO - £2q0- The function 



F(a) - exp 



f fl , da' 

-3 (1+wV)) — 
J\ a' 



(6) 



is determined by the time-dependence of the ratio w{a) between 
the pressure and the density of the dark energy. Equations (01 
and (O show that the observable, i.e. the luminosity distance, is 
a double integral over the dark-energy equation-of-state w(a), 
implying that its dependence on the details of w{a) is quite 
weak. Several articles in the literature have highlighted the pit- 
falls which such a weak dependance produces on the possible 
conclusions that can be drawn on the dar k-energy properties 
dMaor et al.ll2doTl l2002t iBassett et aUl2004l) . 



2.2. Model-independent determination of the expansion 
function 

Instead of specifying a particular Friedmann model and con- 
straining the parameters contained in E{a) as outlined in Eqs. (f5]l 
and (O, our goal is to recover the expansion rate of the universe, 
E(a), as a function of the scale factor a, without assuming any 
specific parameterisation for it. For simplicity of notation, we 
put K — now, allowing us to write /k(w) = w according to 
Eq. (O, which is a first-order approximation even in the case of 
small K + 0. This simplification could be dropped if necessary 
without any change of principle. 

Combining Eqs. (|2), (O and (0), we write the luminosity dis- 
tance as 

c 1 f ' dx cl f'di^ 

H a J a x z E(x) H a J a x l 



x 2 E(x) 

where Ho is the Hubble constant, and we have defined the inverse 
expansion rate e(a) = E^ 1 (a). For the sake of simplicity, we 
shall drop the normalising Hubble length c/Hq in the following 
discussion, thus scaling the luminosity distance by the Hubble 
length. 

Differentiating Eq. (0 with respect to a, and dropping c/Ho, 
we obtain 

(8) 



n ,, . 1 f dx e(a) 

a 1 J a x z a 5 



which can be brought into the generic form of a Vol terra integral 
equation of the second kind for the unknown function e(a), 



?(a) 



, , r dx 

-a D L (a) + A I —e(x) 



(9) 



with the inhomogeneity f(a) = -a 3 D' L (a) and the simple kernel 
K(a, x) - x" 1 . The general parameter A will l ater b e specialised 
to A- a. As detailed e.g. in lArfken & Weberl d 19951) . Eq. © can 
be solved in terms of a Neumann series, 



e(a) = Y J A i e i (a) 



(10) 



where a possible (but not mandatory) choice for the functions e, 

is 



e (a) = f(a) , e„(a) = K(a, f)e n -i (t)dt . 



r- 



(id 



The first guess for eo(a) is equivalent to say that the integral 
or the parameter A in Eq. (0 is small. This crude approxima- 
tion, which is valid in all relevant cosmological cases, is then 
improved iteratively until convergence is achieved. 
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2.3. Application to type-la supernovae 

After application of the empirical relation between duration and 
luminosity, observations of type-la supernovae yield measure- 
ments of the distance moduli p, and redshifts zi for a set of N 
objects, which can be converted into a set of luminosity dis- 
tances D\Sfli) dependent on the scale factors a, = (1 + zd ■ 
For a review about type-la supernovae and their cosmological 
implications see e.g. Leibundgutl (l2001l) . 

As Eq. (0 shows, our analysis requires taking the derivative 
of the luminosity distance with respect to the scale factor. Due 
to measurement errors and scatter of the data about the fidu- 
cial model, it is not feasible to directly differentiate the lumi- 
nosity distance data, since the result would be extremely noisy 
and the estimate of D' L (a) unreliable. Thus, we need to appropri- 
ately smooth the data. We propose to do so by fitting a suitable 
function D^ia) to the measurements D^a,) and approximating 
the derivative in Eq. (0 by the derivative of Dl(o). This choice 
is justified under the assumption that the derivative of the fitted 
data is in fact a good representation of the actual derivative of 
the data. 

For doing this in a model-independent way, it is convenient 
to expand D^{a) into a series of suitably chosen orthonormal 
functions pj(a), 



with 



e ( J\(x)x- 2 dx (18) 

according to Eq. dTTT >. These modes of the inverse expansion 
function can be computed once and for all for any given or- 
thonormal function set \pj(a)}. Due to the linearity of the prob- 
lem, the final solution is then given by 



e(a) = y Cje J \d) . 



(19) 



7=0 



2.4. Error Analysis 

We now show how the errors on the supernova distance mea- 
surements propagate into the expansion coefficients Cj and even- 
tually into the expansion rate. We consider the Fisher matrix of 
the^ 2 function given in Eq. ( fT3l l. 



F tj ee 



2, X 



she 



M 



D L (a) = 2^ c jPj( a ) 



(12) 



dcidcj i 
which is in our case given by 



(20) 



7=0 



The M coefficients c,- in Eq. dT2b are estimated via minimisation 
of the x 1 function 



fc=0 



Pi(a k )pj(a k ) 



crt 



(21) 



X 2 = (Cobs - D(a)f C 1 (D obs - D(a)) , 



(13) 



where D obs is a vector containing the N measured luminosity 
distances, a is a vector of the measured scale factors, and 



where k runs over all supernova measurements and the <x 2 are 
the individual errors on the luminosity distances. By the Cramer- 
Rao inequality, the errors Ac, satisfy 



M 

D(fli) = 2] CjPj(ai) = <Pc)i 

7=0 



(Ac,) 2 > (F- l ) u . 



(22) 



(14) 



is the vector of model luminosity distances to the scale factors a. 

In the final expression of Eq. ( fT4l . P is an N x M matrix M 

with elements P,y = Pj{a,), and c is the M-dimensional vector [Ae(a)l 2 = V 
of expansion coefficients. Assuming that the covariance matrix ^- J 

C~' is symmetric, the set of coefficients minimising^ 2 is 



These errors will propagate into the estimate e(a) of the (inverse) 
expansion function given in Eq. ( fT9] l, 



7=0 



de(a) 



dej 



M 



(Ac/ = J] [e U) («)] (Ac/ . (23) 

7=0 



c = (p T cr l py 1 (p T cr 1 )D 



(15) 



Since the expansion rate is E(a) = l/e(a), its error is finally 
given by 



In this representation of the data, the derivative of the 
luminosity-distance function is simply given by 



[AE(a)Y = 



2 [Ae(fl)] 2 



e 4 (a) 



(24) 



D'Aa) 



Tj C J p 'j (a) 
7=0 



(16) 



thus avoiding the noise which would be introduced by a direct 
differentiation of the data. 

Using the linearity of the integral equation (0, we can now 
solve it for each mode j of the orthonormal function set sepa- 
rately. Inserting the derivative of a single basis function p'Xa) in 
place of D' L (a) in Eq. (0, its contribution to the is given in terms 
of the Neumann series 



„(;) 



(a) = 2" 



< k ef{ 



(17) 



t=o 



3. Application to synthetic and real data sets 

In this section, we demonstrate using synthetic data sets how 
our method performs in two different model cosmologies, an 
Einstein-de Sitter and a standard ACDM model, using simulated 
samples with the characteristics of both current and future sur- 
veys. A brief discussion on the convergence of the Neumann se- 
ries will follow. We also present a toy model with a sudden tran- 
sition in the expansion rate, and show how our method performs 
in this case. Application of the method to the first year SNLS 
data is also presented, as well as the effects of adding high red- 
shift (z > 1) objects. 
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3.1. Illustration: Einstein-de Sitter model 

We start with a simple and unrealistic model cosmology in or- 
der to illustrate the proposed method in detail, i.e. an Einstein-de 
Sitter universe with matter-density parameter Q m o = 1, vanish- 
ing cosmological constant Q.^ = and Hubble constant h - 0.7. 
The expansion function is 



E(d) = a 



-3/2 



e(a) = a 



3/2 



and the luminosity distance is simply 



D L (a) = -(1 - V«) 
a 



(25) 



(26) 



in units of the Hubble radius c/Hq. A suitable choice for the or- 
thonormal function set could start from the linearly independent 
set 



Uj(x) = x 



J/2-1 



(27) 



which can be orthonormalised by the usual Gram-Schmidt pro- 
cedure. The orthonormalisation interval should be [a m i n , 1], 
where a m \ n — (1 + Zmsad is the scale factor of the maximum 
redshift z max in the supernova sample. We thus obtain a set of or- 
thonormal functions \pj(a)}. Projecting the distance in Eq. d26i i 
onto the basis functions, it is straightforward to see that only 
the first two modes po and p\ have non-vanishing coefficients. 
We then use the derivatives of po and p\ to construct the corre- 
sponding Neumann series following Eq. ( TT8b . and from them we 
recover the (inverse) expansion rate. See AppendixlAlfor further 
detail. 




Fig. 1. The reconstructed expansion rate for a simulated sam- 
ple of supernovae in an Einstein-de Sitter universe. The observa- 
tional characteristics of the sample resemble those of the 1st year 
SNLS data. The green shaded area represents the reconstruction 
with 1-cr errors thereof, the blue curve represents the model. The 
bottom plot shows the residuals between the reconstruction and 
the model. 



We applied this method to a synthetic sample of type-la su- 
pernovae in the Einstein-de Sitter universe. The observational 
characteristics of the sample, such as its size, the redshifts and 
the distribution of typical errors of individual measurements, ar e 
adapted to those of the first-year SNLS data (lAstier et alj|2006l) . 
Thus, our sample consists of 120 supernovae up to redshift z — 1 . 



It enables us to determine the expansion coefficients cq and c\ 
with relative errors of order (1-2)%. The reconstructed expan- 
sion rate H(a) - H()E{a) is shown in Figfl] 

The purpose of this simplified example is to show that it is 
possible to achieve a robust and highly accurate reconstruction 
of E(a) when the relevant expansion coefficients can be obtained 
from the data with suitable significance. We now proceed to a 
more realistic application. 

3.2. ACDM model 

We now repeat the preceding analysis with a synthetic sample 
simulated in a standard ACDM universe with Q m0 = 0.3, Qao = 
0.7 and h - 0.7. The expansion function is 



E(a) = (Q t 



i0« 



Q 



AOJ 



1/2 



(28) 



In this case the first two modes of the basis \pj{a)} chosen before 
are insufficient to reproduce D^ia) accurately. When we con- 
sider the true coefficients of the expansion of D^(a), which we 
obtain by projecting it onto the different basis functions, at least 
the first five coefficients are significantly different from zero. 
This is illustrated in Fig. [2] where we compare the model lu- 
minosity distance to its reconstruction using the basis functions 
and the true coefficient of its expansion. If only three or four co- 
efficients are included, the reconstruction deviates significantly 
from the model even at low redshift. 

However, the measurement errors on the data play a cru- 
cial role in this analysis. With current standard data sets, only 
the first three coefficients can be determined significantly, while 
more than just three are needed to achieve an accurate recon- 
struction with the proposed basis functions. We show the re- 
constructed expansion function, obtained including three coef- 
ficients, in Fig. [3] where it is compared to the expansion rate 
of the underlying cosmological model. If we consider only the 
first three modes, all the coefficients are statistically significant, 
although we know from the theoretical model that the recon- 
struction is incomplete. The errors on the first two coefficients 
Co and c\ are of order (1-2)%, increasing to 8% on q. The er- 
rors on higher-order coefficients are larger than the coefficients 
themselves, indicating that they become compatible with zero 
and should therefore be excluded from the reconstruction. 

The precision with which coefficients can be de- 
termined from the data is likely to improve dramati- 
cally with future generation, space-based supernova sur- 
veys such as the Supernova/Accelera tion Probe (SNAP, 
lAldering & the SNAP collaboration! d2004l) ). SNAP is expected 
to measure high-quality light curves and spectra for » 2000 
type-la supernovae in the redshift range 0.1 < z < 1.7. With 
data of such high quality it will become possible to achieve an 
extremely accurate reconstruction of the expansion rate with our 
method. As discussed above, we need at least five coefficients in 
order to reconstruct the expansion rate of an underlying ACDM 
model with the set of functions described above. 

We produced a synthetic data set with SNAP characteris- 
tic s, following the expect ed SNAP redshift distribution reported 
in IShafieloo et al.l d2006l). We also added 25 m ore supernovae 
with z < 0.1, as done in Shafielo o et al.l d2006l) . which are sup- 
posed to be observed by future low-redshift supernova experi- 
ments. Applying our reconstruction technique, we significantly 
constrain the first five coefficients, with errors on the first two 
coefficients being of order 0.1%. The result, obtained using five 
coefficients, is shown in Fig. |4] together with 1-cr errors. 
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Fig. 2. The expansion rate in a ACDM model with Q m o = 0.3 and Qa = 0.7 (solid line) and its reconstruction obtained using the true 
coefficients (dashed line), truncated up to the third (left panel), fourth (central panel) and fifth (right panel) coefficient, respectively. 
The difference between the reconstruction and the model is shown in the bottom panels. When the fifth coefficient is included, the 
two curves nearly coincide up to a — 0.4. 
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Fig. 3. The reconstructed expansion rate for a simulated sam- 
ple of supernovae in a ACDM universe with observational char- 
acteristics resembling those of the 1st year SNLS data. The 
green shaded area represents the reconstruction with 1-cr errors 
thereof, the blue curve represents the model. We made use of 
three coefficients. The bottom plots show the residuals between 
the reconstruction and the model. 



Fig. 4. The reconstructed expansion rate for a simulated sample 
of supernovae in a ACDM universe with the forecast observa- 
tional characteristics of the SNAP experiment. The green shaded 
area represents the reconstruction with 1-cr errors thereof, the 
blue curve represents the model. We made use of five coeffi- 
cients. The bottom plot shows the residuals between the recon- 
struction and the model. 



The choice of the orthonormal function set is in general arbi- 
trary. Obviously, for each underlying model there will be a pre- 
ferred function set, in the sense that the number of coefficients 
required to reproduce the expansion rate is minimal when using 
such a set. It is certainly possible to find a more suitable function 
set for the ACDM model, but since our ultimate goal is to recon- 
struct the expansion rate from the observed data introducing as 
little theoretical prejudice as possible, we are not primarily in- 
terested in finding the most suitable function set to reproduce 
the ACDM expansion rate. Thus, we only made use of the basis 
described in Sect. 13. II 

3.3. Convergence of the Neumann series 

A separate, but related issue is to what power of the parameter A 
we need to follow the Neumann series, or, equivalently, to what 
power k of the scale factor a we need to expand in Eq. ( fTTI i. 
The truncation criterion must again be based on the quality of 
the data. Convergence of the series is achieved at different pow- 
ers k for different redshift intervals. In order to achieve conver- 
gence on the interval 0.5 < a < 1, the series can be truncated 



after k — 4. However, the inclusion of a fourth-order term pro- 
duces a difference to the preceding three orders which is already 
within the error bars, and can therefore be neglected. This trend 
is clearly enhanced when more coefficients are included in the 
reconstruction, since in this case the errors are larger. 

3.4. Recovery of sudden transitions in the expansion rate 

As we have emphasised above, our method can obtain the ex- 
pansion function E(a), or rather its reciprocal e(a), based on a 
representation of the derivative of the measured data. We argue 
that dealing with the derivative of luminosity distance data is not 
expected to cause a major problem, based on the reasonable as- 
sumption that the luminosity distance is a very smooth function. 
As it is evident from Eq. (0, D^(a) is related to the expansion 
function via an integral. Hence, even if E(a) had a peculiar fea- 
ture at some intermediate redshift, this would be smoothed out 
by the integration. 

We address the issue by means of a toy model where the ex- 
pansion function has indeed a sudden transition. We construct 
the toy model starting from the expansion rate of the Einstein-de 



C. Mignone & M. Bartelmann: Measuring the Expansion Rate from Type-la Supernovae 



Sitter model and deforming it by a gentle jump at some interme- 
diate value a* of the scale factor, 



. • arctan [y (a - a,)] + 6 
E(a) = \ a _ 3/2 + j 



(a > a) 
(a < a) 



(29) 



Using Eq. (O, we can obtain the corresponding luminosity 
distance, which is quite smooth and deviates from its Einstein- 
de Sitter counterpart in a way depending on a„. The expansion 
rate and the luminosity distance of this toy model are plotted in 
Fig-El compared to those of the Einstein-de Sitter model. 

Again, we create a synthetic sample of type-la supernovae 
within this model, with the same observational characteristics of 
either SNLS or SNAP, and we apply our reconstruction proce- 
dure. In order to reproduce the transition feature, we need more 
than three coefficients. This is not feasible with SNLS-like data 
because coefficients beyond the third lose significance. With a 
SNAP-like sample instead, the expansion rate can be recovered. 
The results obtained with both synthetic samples, with three co- 
efficients for the SNLS and six for the SNAP case, are shown in 
Fig. [6] together with their 3-cr errors. Figure [6] shows that our 
method can also recover expansion histories with unexpected 
transitions, even though the reconstruction is less accurate than 
that of a perfectly smooth expansion rate. 

We also try to fit this sample to a flat ACDM model and 
explore the parameter space spanned by Q m o and w. The dark- 
energy equation-of-state parameter w is allowed to differ from 
-1, first assuming that it is constant in redshift and then param- 
eterising its time evolution according to 



w(a) = wq + w a (\ - a) . 



(30) 



as proposed bv lChevallier & PolarskH d200ll) and lLinder! d2003l) . 
We find that all the models we considered are capable of pro- 
ducing good fits to the luminosity distance data, but they all 
fail to reproduce the underlying expansion rate when the best 
fit parameters are inserted back into Eq. ©. In most cases, the 
likelihood has more than one maximum, since different combi- 
nations of the considered parameters constrain the two different 
branches of the expansion rate. Unless the time evolution of the 
dark-energy equation-of-state is modelled ad hoc, it is very un- 
likely to reproduce the sudden feature of the toy model in this 
way. However, our method achieves this because the parameters 
involved in our fit trace the relation the between luminosity dis- 
tance and the expansion rate. The different results obtained with 
the usual approach and our method for the fit to the luminosity 
distance and for the expansion rate are displayed in Fig. (0. 

3.5. Application to the first-year SNLS data 

We finally apply our method to the first-year SNLS data 
dAstier et alj 120061) . The sample consists of 71 new super- 
novae observed from the ground with the Canada-France-Hawaii 
Telescope, the farthest being at redshift z = 1.01, plus 44 nearby 
supernovae taken from the literature. Thus, the total sample con- 
tains 115 supernovae in the redshift range 0.015 < z < 1.01. 

As suming a flat, ACDM universe with constant w - -1, 
lAstier et all Ho06) obtained a best fit of Q. mQ = 0.263 ± 0.037. 
Releasing the flatness assumption and adding constraints from 
the baryon aco ustic oscillations (BAO) measured in the SDSS 
dEisenstein et all l2005l) . they obtained Q m0 = 0.271 ± 0.020 
and Qao = 0.751 + 0.082. Furthermore, they investigated mod- 
els with constant equation of state w + -1: assuming flat- 
ness and the BAO constraints, their best-fit parameters are 
Q m0 = 0.271 ± 0.021 and w = -1.023 ±0.087. 




toy model 
Einstein-de Sitter 



0.3 0.4 0.5 0.6 



a 

Fig. 5. The expansion rate of our toy model compared to the 
Einstein-de Sitter one (upper panel), and the corresponding lu- 
minosity distance (lower panel). The parameters for the toy 
model are: a* = 0.7, a = 0.6, y — 11, 8 = 2.3. 



The fit to the luminosity-distance data obtained applying our 
method to this sample, with the orthonormal function set de- 
scribed in Sect. 13.11 is shown in the upper panel of Fig. [8] It 
yields three significant expansion coefficients because the data 
quality especially at high redshift does not allow constraints of 
higher-order modes, as discussed in Sect. 13.21 

The expansion rate reconstructed with our method is com- 
pared in the lower panel of Fig. [8] to that of the best-fit model of 
the SNLS analysis, i.e. a flat ACDM with Q m0 = 0.263. 

3.6. Extending the sample beyond z=1 

Another interesting problem concerns what could improve the 
performance of the method. Clearly, both a larger sample of su- 
pernovae and a better accuracy in the individual measurements 
would help reducing the errors on the coefficients, which would 
eventually enable a significant estimate of more coefficients, and 
thus a more precise reconstruction of the expansion rate. From 
a mathematical point of view, adding more objects and reducing 
the uncertainties are equivalent: a sample four times larger than 
another yields the same results obtained with the smaller sample, 
if its error bars were to be reduced by one half. However, since 
the measurement accuracy cannot be indefinitely shrunk below a 
given limit, because of systematic uncertainties, a long run strat- 
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Fig. 6. The expansion rate of our toy model (blue curve) and 
its reconstruction, with 3-<x errors thereof (green shaded area), 
obtained from a SNLS-like data set with three coefficients (up- 
per panel), and from a SNAP-like data set with six coefficients 
(lower panel). The bottom plots show the residuals between the 
reconstruction and the model. 



egy to make best use of the method would be to increase the size 
of the sample. 

Extending the sample to higher redshift can also help reduc- 
ing the errors on the estimated coefficients. We investigated the 
issue by means of an extremely simplified example: we added to 
our previously described simulated ACDM SNLS-like sample 
of supernovae up to z - 1 an additional set of 20 objects, uni- 
formly distributed between z - 1 and z - 1.7. The excess in the 
Fisher matrix due to the inclusion of the higher-z sample influ- 
ences the errors on the coefficients, reducing the first by ~ 10% 
and the following ones by ~ 20%, whereas if we only added 20 
more objects with z < 1 the gain would be smaller. However, 
this still does not allow more than three coefficients to be signifi- 
cantly pinned down, with the orthonormal function set described 
inSect.|3~T1 

We also apply our method to a supernova sa mple which 
extend s beyond z — 1, namely the one compiled by iDavis et al.l 
d2007l) . inc luding the comb i ned ESSENCE/SNLS/nearby 
dataset from Wood-Vasev et al.l (|2007) and the HST data from 
iRiess et all d2007l) . It contains 192 supernovae, of which 15 
with 1 < z < 1.75. Although this sample contains more objects 
than the SNLS and extends to higher redshifts, the quality of 
the reconstruction achieved is not better than the one obtained 



Fig. 7. Upper panel: The luminosity distance of the toy model 
(cyan curve) together with the SNAP-like simulated sample 
(black points), compared to our fit (blue dashed curve) and three 
other cosmological fits (red curves). The bottom plot shows the 
residuals between the different fits and the model. Lower panel: 
The expansion rate of the toy model, of our reconstruction and 
of the other models. The different red curves correspond to three 
different fits to a flat ACDM: in case #1 (red dashed curve) we 
impose w = — 1 and let only Q m o vary, in case #2 (red dotted 
curve) we let both Q m o and w (constant in redshift) vary, and in 
model #3 we let Q m o, wo and w a vary, according to Eq. d30l ). 

using the SNLS data set. In fact, the errors on the coefficients 
are slightly larger, because the individual uncertainties on 
the distance moduli in the extended ESSENCE sample are 
significantly higher than the SNLS ones (at least for z < 0.8), 
due to the different way luminosity distances are estimated 
from the photometric data by the two groups (a r eview on the 
different supernova light-curve fitters is given in IConlev et al.l 
d2007h V 

4. Discussion 

We are proposing a method to constrain the expansion func- 
tion of the universe without assuming any specific model for 
Friedmann-type expansion. If the universe is isotropic, homo- 
geneous and simply connected, it is described by a Robertson- 
Walker metric. Cosmological measurements can generally only 
constrain one of two functions of time, the expansion rate and 
the growth of structures. We have shown here how the expan- 
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Fig. 8. Upper panel: The 1st year SNLS sample and our fit for 
the luminosity distance. Lower panel: Our reconstruction of the 
expansion rate and 1-cr errors thereof (gray shaded area), com- 
pared to the ex pansion rate of a A CDM model with the best fit 
parameter from lAstier et all d2006h . Q m0 = 0.263 ± 0.037, and 
1-cr errors thereof (red curves). 



The drawback of the transformation to a Volterra integral 
equation is that the derivative of the luminosity distance with 
respect to the scale factor is needed to start the Neumann series. 
Derivatives of data are notoriously noisy and should be avoided. 
We propose to expand the luminosity distance into an initially 
arbitrary orthonormal function set, fit its expansion coefficients 
to the data and then use the derivative of the series expansion 
instead of the derivative of the data. Suitable orthornormal func- 
tion sets can be constructed by Gram-Schmidt orthonormalisa- 
tion from any linearly independent function set. The only con- 
dition so far is that the number of coefficients required to fit the 
data should be minimal. 

Once the orthonormal function set is specified, the Neumann 
series can be constructed beforehand for all its members. The 
measured coefficients of the series expansion directly translate 
to the solution for the expansion function. The convergence cri- 
terion for the Neumann series is determined by the data quality, 
as is the number of orthonormal modes in the series expansion 
of the data. 

Applications to synthetic data samples of increasing com- 
plexity are very promising. In particular, we showed that an ex- 
pansion function containing a sudden transition can be faithfully 
recovered by our method provided sufficient quality of the in- 
put data. We also apply our method to the first-year SNLS data 
and show the expansion function recovered, and notice how the 
inclusion of higher-redshift (z > 1) objects could improve the 
performance of the method. 

In future studies, we shall investigate how our method can 
be used to determine the single expansion function underlying 
all cosmological measurements used. 
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sion function can be observationally constrained without refer- 
ence to specific assumptions on the time evolution of the terms in 
Friedmann's equation and their parameterisation in terms of den- 
sity parameters. The issue may become important in search of 
constraints for a dynamical dark energy component, whose be- 
haviour is so far only very poorly known. Since it is unclear how 
its energy density contribution to the cosmic fluid may change 
in time, any guessed parameterisation may be erroneous, hence 
a parameter-free recovery of the cosmic expansion rate may turn 
out advantageous. 

We demonstrate our method here using the luminosity- 
distance measurements obtained from type-la supernovae as a 
model. Since the luminosity distance is a cosmological observ- 
able depending on space-time geometry only, the dynamics of 
structure growth does not enter yet. The method proceeds in two 
essential steps. First, the integral relation between the expan- 
sion function and the luminosity distance is transformed into a 
Volterra integral equation of the second kind. Under the relevant 
conditions, its solutions are known to exist and to be uniquely 
described in terms of a convergent Neumann series. In other 
words, the method is guaranteed to return the unique expansion 
rate of the universe within the accuracy limits allowed by the 
data. 



Appendix A: Exact Solution for the Einstein-de 
Sitter case 

Here we show how to construct the (inverse) expansion rate of the Einstein- 
de Sitter model from the first two modes of the function set obtained applying 
Gram-Schmidt orthonormalisation to the set of linearly independent functions 
specified by Eq. (27). The first two modes are 



1 1 

Po(x) = — — - 



pi(*) = -= -= + — 
Vc\ V* * 



(A.l) 
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It is straightforward to see, by projecting the distance in Eq. )26t onto the basis 
functions, that only the first two modes are needed, i.e. 



Ol(«) = ^ CjPj(a) , 



(A3) 
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where Cj = J D^(a)pj(a)da stands for the 7-th true coefficient of the expan- Daly, R. A. & Djorgovski, S. G. 2003, ApJ, 597, 9 



sion. In this case £q = 2(1 + 2/3) -\/a and f 1 
From the derivative of po(a), 
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we can construct the zero-th order Neumann series following Eq. 1111 : 
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Up to second order, the zero-th order Neumann series for the (inverse) expansion 
rate is 
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Again, from the derivative of p\(a) 



P ' W = -4(^ + ?) 



we can construct the first-order Neumann series: 
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The first-order Neumann series up to second order thus reads 
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Now we can employ Eqs. JA.6I and jA.lOt and the true coefficients of the 
expansion, and recalling the relations in Eq. IA.2K we recover the inverse expan- 
sion rate for an Einstein-de Sitter universe: 
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